Equilibrium points and their stability of COVID-19 in US

This study aims to develop an advanced mathematic model and investigate when and how will the COVID-19 in the US be evolved to endemic. We employed a nonlinear ordinary differential equations-based model to simulate COVID-19 transmission dynamics, factoring in vaccination efforts. Multi-stability analysis was performed on daily new infection data from January 12, 2021 to December 12, 2022 across 50 states in the US. Key indices such as eigenvalues and the basic reproduction number were utilized to evaluate stability and investigate how the pandemic COVD-19 will evolve to endemic in the US. The transmissional, recovery, vaccination rates, vaccination effectiveness, eigenvalues and reproduction numbers (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${R}_{0}$$\end{document}R0 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${R}_{0}^{end}$$\end{document}R0end) in the endemic equilibrium point were estimated. The stability attractor regions for these parameters were identified and ranked. Our multi-stability analysis revealed that while the endemic equilibrium points in the 50 states remain unstable, there is a significant trend towards stable endemicity in the US. The study's stability analysis, coupled with observed epidemiological waves in the US, suggested that the COVID-19 pandemic may not conclude with the virus's eradication. Nevertheless, the virus is gradually becoming endemic. Effectively strategizing vaccine distribution is pivotal for this transition.

www.nature.com/scientificreports/Susceptible-Infected-Recovered (SIR) model, a staple in prior research, has been instrumental in simulating COVID-19 spread within communities [12][13][14] .Researchers have expanded the model to include an 'Exposed' (E) compartment, thereby evolving it into the SEIR model for a more accurate representation of the disease's progression [15][16][17][18][19] .However, since December 2020, the USFDA's emergency authorization of the Pfizer and Moderna mRNA vaccines marked a turning point 20 .With Phase 3 clinical trials indicating approximately 95% efficacy for both vaccines in preventing COVID-19 21 , governments and communities have been actively promoting vaccination.By the end of 2021, around 73% of the U.S. population had received at least one vaccine dose, with 62% fully vaccinated 22 , profoundly affecting the virus's transmission dynamics.This significant shift highlights the urgency to incorporate vaccination data into epidemic models 23,24 .On the other hand, the traditional mathematical models rely on a set of differential equations, merely estimating parameters from real data through numerical solutions of these equations falls short in fully capturing the dynamic nature of COVID-19 transmission and in formulating effective strategies to combat the virus.Stability analysis, using real-world data and often neglected in COVID-19 epidemic models, is crucial for understanding the virus's dynamic behavior.COVID-19's tendency to oscillate between stationary and non-stationary states leads to multiple epidemic waves, a complex pattern that demands thorough analysis.
To tackle these challenges, we propose an augmentation to the SEIR model by integrating a vaccination component, resulting in the SEIRV model.Our primary objective involves conducting comprehensive stability analyses of the SEIRV model within the context of COVID-19 10,11,22,[25][26][27] .We intend to delineate both the diseasefree and epidemic critical equilibrium points and subsequently employ advanced next-generation matrix methods to compute the basic reproduction numbers pertaining to these critical points: R 0 for the disease-free state and R end 0 for the epidemic state.Employing stability theory, we aim to derive conditions conducive to stable states within the SEIRV system.Moreover, our endeavor includes the identification of attractors defining the spectrum of parameters governing stable states within the dynamic system of COVID-19 28 .Leveraging this multifaceted stability analysis, our overarching aim is to provide critical insights to aid in devising strategies aimed at mitigating and ultimately arresting the transmission of COVID-19 in the United States.

Data sources
The number of cases and deaths for each state from January 21, 2020 to December 12, 2022 was downloaded from https:// github.com/ nytim es/ covid-19-data.The vaccination data, including the number of vaccines distributed, the number of people who received at least one shot of vaccinations and full vaccination for each state from January 12, 2021 to December 12, 2022 were downloaded from https:// ourwo rldin data.org/ us-states-vacci natio ns.

SEIRV compartment model for COVID-19 transmission dynamics
We constructed a compartmental model based on a deterministic system of nonlinear differential equations for COVID-19 transmission dynamics, which further took into account vaccination.The dynamical system was established based on the assumptions: 1. We assumed that the total population could be classified into several compartments, and the transition of individuals from one compartment to another depends on the stage of the disease; 2. We assumed that if susceptible individuals exposed to an infectious source, they are not infectious until they have symptoms after a latent period; 3. We assumed that the infected individuals will recover or die after an infectious period, however, for model simplicity we put them into one group which is called recovered; 4. We assumed a constant recruiting rate and natural death rate.
Based on these assumptions, we established the model including the susceptible ( S(t) ), exposed ( E(t) ), infected ( I(t) ), recovered ( R(t) ), and vaccinated ( V (t) ) individuals, each one forming a compartment.The model employed a constant recruiting rate (births) to the susceptible individuals ( S(t) ), natural death rate µ , transmission rate β, vaccination rate α , incubation rate γ , the probability of the recovery or death δ , and vaccine inefficiency σ .The detailed definition of these populations and description of the transmission relationship among them have been concluded in Table 1.
We plotted Fig. 1 to depict the transitional relationships among the various states within the model.Based on the above assumptions and concepts, the rates of change of the five populations were governed by the following system of nonlinear ordinary differential equations, which constituted the SEIRV model used in our study: www.nature.com/scientificreports/Summarizing the above equations, we obtained From Eq. ( 6), we can rewrite N(t) as We assumed that N(0) ≤ � µ , so we obtain N(t) is bounded by � µ for all t, Thus, we showed that all variables S(t), E(t), I(t), R(t), V (t) were bounded in the region , and more details could be found in Supplementary Note S2:

Reproduction number
The basic reproduction number, denoted as R 0 , is an important threshold quantity which determines whether the COVID-19 will continuously spread in the population or disappear.It is defined as the average number of secondary cases produced by one infected individual introduced into a population of susceptible individuals 29 .
For the set containing all infected individuals ( E(t) and I(t) ), we defined There is

Let
The basic reproduction number R 0 was defined as the spectral radius of the next-generation matrix F ′ (X)V ′ −1 (X) as the following, and more detailed calculation was illustrated in Supplementary Note S3.1 29 : The reproduction number R 0 is used to measure the transmission potential of a disease.Intuitively, we can expect that if R 0 < 1 then the number of new cases of COVID-19 will decrease, and the number of new cases will increase if R 0 > 1.

Steady sate analysis of SEIRV for COVID-19 transmission dynamic system
Predictions of multistable structural dynamics are paramount to the development and implementation of vaccination and population public health intervention for controlling the spread of COVID-19 under highly uncertain biological, economic, political, and environmental perturbation.Although a direct numerical solution to differential equations can be used to calculate the steady state, analytic analysis of the steady state can reveal how the parameters affect the steady state and enable prediction of near-and far-from equilibrium response.Such analyses are invaluable for understanding the endemic waves of COVID-19 and developing strategies to mitigate or even end the pandemic.
COVID-19 transmission dynamic system is an autonomous system in which the independent time variable t does not explicitly appear in the differential equations.We will focus on steady state and stability analysis of autonomous systems to investigate the qualitative dynamic behavior of COVID-19 dynamic systems.The central issue in stability analysis is to identify isolated critical (equilibrium) points of COVID-19 transmission dynamics.The isolated critical points can be classified as disease-free (number of new cases is zero) critical point and endemic (new cases are not zero) critical point.The two classes of critical points are given as follows (For details, please see Supplementary Note S3.2).

Classification of critical points
Stability analysis for the general nonlinear dynamic systems is complicated.In this paper, we will focus on isolated critical point and almost linear systems.A system is called almost linear at a critical point if the Jacobian matrix of linearized system at an isolated critical point is invertible.We assumed Jacobian matrix of the nonlinear dynamic system (1-5) is invertible, and could be derived as where Its characteristic polynomial is where are the corresponding eigenvalues.Using Routh-Hurwitz stability criterion, we could obtain the positivity or negativity of the eigenvalues, which could help us study the properties of the critical points which determine whether the system is stable or unstable or whether the number of new cases decreases or increases.Disease-free and endemic critical points are separately investigated and briefly presented (for details, please see Supplementary Note S3.3).
. Then, in summary, the endemic equilibrium point can be classified into three cases: (1) when R end 0 < 1 , the endemic equilibrium point is classified as an asymptotically stable node; (2) when R end 0 = 1 , the endemic equilibrium point is classified as an unstable node; (3) when R end 0 > 1 , the endemic equilibrium point is classified as an unstable saddle point.

Parameter estimation
We can use the steady states of the COVID-19 transmission dynamic system to estimate the parameters in the nonlinear differential equation models and then use the estimated parameters to identify and classify the equilibrium points.Let ( 16)

Ethical approval
This study did not involve individual information, so there was no requirement for informed consent.

Estimation of parameters
The data we employed started from January 12, 2021, and ended at December 12, 2022.To examine endemic equilibrium points, we identified three distinct steady state periods: (1) April-July, 2021; (2) March-May, 2022; and (3) September-November, 2022 (Fig. S1).The daily reported new infection cases were utilized to calibrate the COVID-19 model, enabling precise estimation of the model parameters.The parameters are estimated by minimizing the Eq. ( 21).We applied the black-hole optimizer (BHO) approach, which is considered as a recent metaheuristic optimization technique and is used to solve continuous optimization tasks by adapting physical behaviors in evolution rules [30][31][32][33] .For all the states, we assumed a same and constant birth rate N and natural death rate µ .The other five parameters α, β, γ , δ, and σ are derived for each state individually using their respective observed number of new cases and total population size.In order to avoid the impact of random initial setting on the results, for each state, we performed estimation of 100 times and selected the parameters that minimized the loss of Eq. ( 21).Within each estimation process, the maximum number of iterations was set as 5000.Figure 2 displayed the estimated median values and standard deviations of the parameters α, β, γ , δ, and σ across three steady periods for all 50 states in the US.We observed that comparing these three steady state periods, the trends for parameters β, σ , γ , and δ showed a decreasing pattern, with only parameter α increasing.For all parameters, the changes from the time periods April-July, 2021 to September-November, 2022 were significant using paired wilcoxon signed-rank test ( β : p < 0.0001 , σ : p = 0.0021 , γ : p = 0.0008 , δ : p < 0.0001 , α : p = 0.0003 ).The decreases in both transmission rate β and incubation rate γ indicated a significant shift in the dynamics of the disease spread, suggesting that the disease was spreading more slowly, likely due to effective public health interventions.Meanwhile, the decrease in vaccine inefficacy rate σ coupled with an increase in vaccination rate α reflected a positive development in vaccine-related defenses.Vaccines became more effective and more widely administered, boosting population immunity.In summary, the changes observed in the parameters suggested a gradual reduction in COVID-19 spread rates in the US, likely attributable to the escalating vaccination efforts.

Eigenvalues, reproduction number at the endemic equilibrium points and stability analysis
To evaluate the stability of COVID-19's dynamics across the 50 states in the US, we analyzed the eigenvalues of the Jacobian matrix at endemic equilibrium points.Eigenvalues could help determine the stability of the virus' spread, assessing whether the virus's spread will be killed or persist within a population.Our system had a total of five eigenvalues, with theoretical stability analysis revealing that one was consistently negative.Notably, in all three steady periods, our calculations revealed that at least one eigenvalue was positive in each state (Table S1), indicating that the COVID-19 dynamics across the US have not yet stabilized.However, there was a notable decrease in the largest eigenvalue in 32 states from the first (April-July, 2021) to the third steady period (September-November, 2022), while 3 states exhibited no significant changes during these periods.The reproduction number R end 0 at the endemic equilibrium point indicates the average number of people that one infected individual will infect in a fully susceptible population, and can be used to characterize how far the dynamic system is from the stable state.Values of R end 0 exceeding 1 suggest a potential epidemic, while values below 1 indicate a possible decline in the disease; if they stabilize around 1, COVID-19 may become endemic, maintaining predictable levels within the population.Table S2 presented the reproduction number R end 0 for all three steady periods across 50 states in the US, revealing that the number of states with R end 0 less than 1 were 2, 4, and 19 in the first, second and third steady period, respectively.Figure 3 featured a map showing the distribution of R end 0 values across these periods in the US.Notably, from the result of the third steady period, most states in the Middle East were expected to reach a steady state in the latest COVID-19 wave in the US.A comparison between Tables S1 and S2 indicated that 16 states with an R end 0 below 1 in the third steady period corresponded to those showing a decrease in the largest eigenvalues from the second to the third period.( 20)

Impact of the parameters on the eigenvalues
To illustrate the impact of the parameters β, α, σ , γ , and δ on the largest eigenvalue of the COVID-19 dynamics, we plotted Fig. 4 to demonstrate the relationship between these parameters and the largest eigenvalue in the COVID-19 dynamic systems.Our observations revealed a positive correlation of the largest eigenvalue with β , σ , γ in all three steady periods, and with δ in the second and third periods.In contrast, there was a negative correlation with the parameter α .Therefore, to lower the largest eigenvalue, and thereby the disease's spread, it's necessary to decrease the rates of transmission β , reduce the time individuals that spent in infectious state γ , increase the recovery rate δ , increase the vaccination rate α , and improve vaccine effectiveness 1 − σ .This combination of strategies would collectively reduce the spread of the COVID-19 within the population.

Approximate attractor analysis and impact of the parameters on stability
An attractor is defined as a set of states toward which a system attempts to return after perturbation 34 .Rather than examining the attractors of the state variables I, V , S, E , and R in our system, we focused on identifying the parameter regions that govern the stability of COVID-19 dynamics and the effects of parameter variations on this stability.Our model incorporated five estimated parameters, presenting a challenge in simultaneously depicting and visualizing the stability region encompassing all these parameters.Since the vaccination rate α and vaccine inefficiency rate σ jointly contribute to the transmission dynamics of COVID-19, we investigated the impact of parameters α and σ together on the stability of COVID-19 dynam- ics. Figure S2 showed the impact of simultaneously varying parameters α and σ while maintaining the current values of other parameters unchanged on the stability of COVID-19 dynamics in the third steady period of September-November 2022 across 50 states in the US.In this figure, the stability region of the parameters was represented in blue color, while the instability region was shown in red color.An upward movement along the vertical axis represented an increase in the vaccination rate α , and a leftward movement on the horizontal axis signified a decrease in the vaccine inefficiency rate σ or, conversely, an increase in vaccine efficacy.Our observa- tions from Fig. S2 indicated that Arkansas (AR), Kentucky (KY), Idaho (ID), and Georgia (GA) had the largest stability regions for α and σ , whereas Alaska (AK), Florida (FL), Louisiana (LA), Connecticut (CT), and Texas (TX) had the smallest.The above finding offered crucial insights that enhancing vaccination rates can compensate for low vaccine efficiency, and similarly, improving vaccine efficiency can balance a lower vaccination rate to maintain stability.
Subsequently, we examined the combined effect of parameters γ (incubation rate) and δ (probability of recov- ery or death) on the stability of COVID-19 dynamics.Figure S3 displayed the impact of simultaneously changing γ and δ , while holding the values of other parameters constant, on the stability of COVID-19 dynamics during the third steady period.Our analysis indicated that increasing δ and decreasing γ would nudge the COVID-19 dynamic system closer to stability.In a similar vein, Figs.S4 and S5 demonstrated the consequences of changing parameters β and γ , and β and δ , respectively, under the same conditions, on the stability of COVID-19 dynamics across the states.Furthermore, Fig. S6 illustrated the effects of altering β, γ , and δ simultaneously, while keeping other parameters constant, on the stability regions in selected states.It specifically represented states with the www.nature.com/scientificreports/largest and second largest stability regions, Michigan (MI) and Kansas (KS), against those with the smallest and second smallest stability regions, Nebraska (NE) and Washington (WA).In summary, these figures showed that implementing strategies decreasing transmission rate β , decreasing incubation rate γ , and increasing the prob- ability of recovery or death δ simultaneously, or any two of them, would drive the dynamic system of COVID-19 toward stability region.Biologically, decreasing β meant fewer individuals are infected over time, lowering the spread of the virus.Reducing γ would extend the time before exposed individuals can infect others, reduce the rate of infection, and provide a chance to contain the virus through early intervention, such as isolation or treatment.Increasing δ would lead to quicker recuperation, shrinking the pool of infectious When these changes occur together, or even just two of them, it would slow the disease transmission and guide COVID-19 dynamics toward stability.

Conclusions
At present, the most pressing and widespread question concerns the end of the COVID-19 pandemic 35 .While some scientists have optimistically forecasted the pandemic's conclusion in 2022 36 , others caution that it is far from over [37][38][39] .In response to this debate, our paper focused on developing a method to analyze the stability of COVID-19 transmission dynamics.This approach serves as a tool to evaluate whether the pandemic is transitioning into an endemic phase.
There is no universally agreed-upon definition of a pandemic 40,41 .Pandemic and endemic are terms rooted in systems dynamics.In our study, we interpreted "endemic" as the "stable states" of the COVID-19 dynamic system.Generally, protection against COVID-19 infection relies on vaccination and natural immunity, which of both are important to be considered in modelling, but the data for natural immunity is missing.In this paper, we proposed an innovative approach, which embodied in the integration of vaccination into the SEIR model to study the stability of the COVID-19 dynamics.By accounting for vaccinated individuals as a separate compartment, we could capture the nuanced impacts of vaccination campaigns on disease dynamics, offering insights into anticipate the effects of vaccine rollouts on the transmission of the virus.This advancement not only enhances the model's ability to simulate real-world scenarios, but also provides a valuable tool for public health planners to optimize vaccination programs and mitigate the spread of the disease.Specifically, we formulated an objective function to ensure the steady states of dynamics which was then applied in estimating the parameters driving the current transmission dynamics of COVID-19.We derived eigenequations to determine eigenvalues, as well as reproduction number R 0 at disease-free equilibrium point and R end 0 at the endemic equilibrium point, to assess the stability of COVID-19 dynamics under these conditions.Our analysis led to two significant conclusions.The first was that all 50 states in the US were in unstable states during all three steady periods.The second conclusion, however, indicated that the COVID-19 dynamics are moving towards stability.Analysis of real data revealed that the three transition-related parameters β, γ and δ decreased while the vaccination rate α and vaccine efficiency 1 − σ increased, progressing from the first steady period through the second one to the third period.Again, the results also showed that both the largest eigenvalue of the stability analysis matrix and the reproduction number R end 0 in three steady periods decreased consistently across the three steady periods.Both vaccination and nonpharmaceutical interventions (NPIs) are critical in curbing the spread of COVID-19.Theoretically, by exploring the parameter space of vaccination-related factors α, σ and NPI-related factors β, γ , δ , we can define the stability region of these parameters during steady state periods.In this paper, we have calculated the stability regions for α, σ , and β, γ , δ , offering insights into transitioning COVID-19 from a pandemic to an endemic state.Enhancing vaccine efficacy, such as through the use of mucosal vaccines, coupled with an increased vaccination rate and the implementation of appropriate NPIs, can facilitate this change from pandemic to endemic.
In conclusion, our research established a novel SEIRV model and thus provided compelling evidence that COVID-19 will not conclude with the virus's elimination but is transitioning towards an endemic stage in the United States.The stability analysis, aligned with observed epidemic patterns, underscored this inevitable shift.Critically, our findings emphasized the importance of strategically designed vaccine deployment as a determinant in this evolution.Effective vaccination strategies are pivotal not only in mitigating the immediate impact of the pandemic but also in shaping the long-term coexistence with the virus.This research contributed valuable insights into the trajectory of COVID-19 and served as a guide for policymakers to formulate responses that could lead to a manageable endemic state, ensuring public health resilience in the face of this ongoing global health challenge.
While our model represented an improvement over previous research, there remained room for further refinement.For example, the rapid mutation rate of viruses, leading to the emergence of new virus variants and reinfection, demanding the model more accurately mirror the real-time complexities of the pandemic.Future research should explicitly incorporate the emergence of new variants with potentially altered transmission characteristics and immune escape capabilities into the model.Additionally, the model needs adapt to the fact that recovered individuals can become susceptible again.Such improvements would enable the model to provide more precise simulations and predictions of transition dynamics of COVID-19, which are crucial for crafting effective and timely public health responses to ongoing and future outbreaks.
The evolution of the virus remains unpredictable.The transition dynamics of COVID-19, influenced by factors such as virus variant evolution, the emergence of new variants, vaccination, and the unknown aspects of natural immunity, are complex and uncertain.There is no straightforward path from a pandemic to an endemic state.As editorial of "nature" pointed out "there won't be a single 'exit' wave to mark the lifting of pandemic restrictions.Further waves of infection and death are likely to follow, either from new variants that arise in the population, or from variants imported as the country opens its borders to visitors" 41 .The US has already experienced five major COVID-19 waves and may face more.However, the encouraging news is that COVID-19 in the US is progressively approaching an endemic stage.

Figure 2 .
Figure 2. Box plot for the parameters β, α, σ , γ , δ in three steady periods (April-July, 2021; March-May 2022; September-November, 2022) across 50 states in the US.(a) Box plot for the transmission rate β , representing the rate at susceptible compartment transmit to exposed one, (b) box plot for the vaccination rate α , representing the rate at susceptible individuals become vaccinated, (c) box plot for the vaccine inefficiency σ , representing the rate at vaccinated people become exposed, (d) box plot for the incubation rate γ , representing the rate at exposed population moves to infected one and (e) box plot for the probability of the recovery or death δ , representing the rate at infected persons become recovered.

Figure 3 .
Figure 3. Map of distributions of the reproduction number R end 0 in three steady periods across 50 states in the US.This figure was generated using ArcGIS 10.6 (ESRI, Redlands, California USA; URL: http:// www.esri.com/).

Figure 4 .
Figure 4.The relationship between the parameters and the largest eigenvalue of the COVID-19 dynamic systems.(a) The relationship between transmission rate β and the largest eigenvalue, (b) the relationship between vaccination rate α and the largest eigenvalue, (c) the relationship between vaccine inefficiency σ and the largest eigenvalue, (d) the relationship between incubation rate γ and the largest eigenvalue, and (e) the relationship between probability of the recovery or death δ and the largest eigenvalue.

Table 1 .
Description of notations and parameters in the model, the first column denotes the symbols and the second one denotes the description.The probability of the recovery or death due to COVID-19, the rate at infected persons become recovered or dead caused by COVID-19 σ The vaccine inefficiency, the rate at vaccinated people become exposed, 1 − σ represents the population vaccine efficacy

Figure 1. Flow
diagram of the dynamic transmission model of COVID-19.